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Variation of Parameters Methods 
with Application to BGI4 

Robert G. Gottlieb 

McDonnell Douglas Space Systems Company 
Engi.eering Services Division 


SUMMARY 

There are a number of variation of parameter methods found in the literature that use ficti- 
tious tune, s, defined by the Sundmann tra ns for m ation r * dt/ds as independent variable. 
The most notable being KS, Kustaanhdmo and SliefeJ [1]J2], Burdet [3], BG14, Band and 
Gottlieb [4], and BG14& Bond[12]. There are two problems that wise when fictitious tune 
is used as independent variable. The first is dm stopping on time becomes an iterative pro- 
cess, and the second, more serious, problem is that, without care, the methods mentioned 
above may cease to work altogether for certain classes of problems. In this papa- a simple 
derivation of several variation of parameters methods using universal variables is given. A 
schematic i nt e rpr et ation of variation of parameters is given which points out the r*nse of 
the failure of previous methods. It is then shown how the independent variable (fictitious 
time) and other constants may be chosen so that stability and accuracy may be maintained 
even for those classes of problems that caused previous methods to fail. Augmenting the 
energy variable all the way to the Jacobi constant as was done in [4] is postponed until the 
end of the paper because it was felt dm it would distract the reader from the issue of 
interest, namely, the optimal use of fictitious time. In addition, it is shown how to handle 
multiple vehicles and constant time steps. The use of constant time steps in conjunction 
with BG14 was found to be as accurate as fictitious time alone for near circular orbits. 
Non circular orbits were not tested. 


INTRODUCTION 

This work was prompted by the discovery of two classes of problems that caused the varia- 
tion of parameters methods KS. and BG14 to fail, in the sense that the computer programs 
would either stop producing output, (V would produce erroneous output. The two classes of 
problems were long term low thrust to escape, and reentry to the ground. It was found 
that restarting these problems at intermediate points caused them to proceed to a solution. 
This indicated that the difficulty in gating them to work without restarting was in some 
sense related to either the size of the increasing independent variable (secular terms), or the 
fixing of the origin, or both. Bond [12] suspected that the demination of secular terms in 
the original formulation [4] would effect a cure and proposed BG146. Although it is stated 
in [12] that the removal of secular terms cures the instability problem, further experimenta- 
tion has shown that, although they run a little further, the long term low thrust to escape, 
and reentry to the ground problems cause BG145 to fail just like the others. This meant 
that restarting was the answer, but a total restart would mean a loss of accuracy. This is 
true because redundant variables are carried, and a total restart would mean that some of 
them would have to be ignored. Therefore, this paper examines how to restart in such a 
way as to lose no information. That is to say full use is made of the redundant variables, 
and none is recomputed. The way to do this is motivated by examining what is going on in 
the derivation of a variation of parameter method. This restart without loss of information 
is what is referred to as the optimal use of fictitious time. It turns out that this feature 
alone, without concent for secular terms, cures the numerical problems, mentioned above, 
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that plagued previous implementations of BG14 and by inference KS. 
DERIVATION OF EQUATIONS 


The differential equation for perturbed two-body motion is 

= _^) +? 
r J ar 


0 ) 


where Y is the position vector of one of the bodies with respect to the other, p is the gravi- 
tational parameter, f is the vector of all perturbing accelerations. V(Yj) is the perturbing 
potential, and 4 is the vector of all perturbing accelerations not derivable from a potential. 
Although eq(l) is easily written down, and is seen often, it is difficult to integrate accu- 
rately for long periods. The problem is the term. This term causes eq(l) to be non- 

r 3 

linear, and to have a singularity at r = 0. These objections can be removed in two steps. 
The first step is to make a change in independent variable from time, t. to fictitious time, s. 
where s is defined by 


This is known as the Sundmann transformation, ami was first used in a theoretical investi- 
gation of die three-body problem [5]. Using the chain rule for differentation, aid denoting 
differentation of ( ) with respect to s by ( )', results in 


also, since r 2 = Y - Y 


Y’ -Y t‘ = Yr 
t' = r l' = rr 


r' = rr =rr 


( 3 ) 

( 4 ) 

( 5 ) 


From eq(3) and eq(4) 

?" = r¥+rfF=r¥+(r , F)F (6) 

The second step is to replace the (YY)Y term in eq(6). It was evidently Sperling [6], 
who first removed the (Y Y)Y term by using the Laplace or perigee vector, p , of two- 
body motion, where 




(Y Y)Y - (Y ■ Y)Y Y 


Therefore, using eq(I) and eq(7) in eq(6) results in 
Y" = r 2 J * - \ip + 


Y Y- 2a 


( 7 ) 

( 8 ) 


The coefficient of Y is easily recognized as twice the energy integral of Keplerian motion, 
so, defining 


a* sY Y 


_2a 

r 


( 9 ) 


eq(8) may be written 

Y" = a t r+r 2 f- V ip (10) 

The introduction of the Laplace vector linearizes the Y " equation when f - 0. The real 
advantage of Sperling’s insight will be apparent later. 


Now, assuming X and V are arbitrary constants, add Xf* and V to eq(10), and subtract XY 
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and V from eq(10). This will leave t " unchanged no matter what X or V are. These con- 
stants are introduced at this point in anticipation of being able to choose them such that the 
integrability of the variation of parameters equations (developed in the next section) is 
improved. Re-writing eq(10) 

r w = a 4 F*+» J /-pjr + Ar-Xy»+tf-tf (li) 

Now, defining 
o«aj + X 

r-ujf+tf 02) 

@ m r 2 /* - XT* 

f mQ + V 
eq(U) becomes 

r" = Qj>+f -t (13) 

Using eq(5) 

r"=(f > i* + r f)r (14) 

Now, using eq(l), eq(9). and the definitions of a and @ in eq(14) results in 

r" = or + p + r £? (15) 

where r s — . At this point it is informative to examine the derivatives of a and t. Start- 
r 

ing with a, and recalling that X is a constant 

a! = ra = 2jrf*t+2^r (16) 

Using eq(l) and eq(S), a' very easily reduces to 

a! = 2(7*' f ) (17) 

The derivative of £ is 

• » 

t’ = ri=r2tf*i^+(?-7*)i*-(7*i*)r-(7*-¥)7*- ( 7 * ?)? - \i- + \iXrr ( 18 ) 

r r J 

Simplifying and using eq(l) and eq(S), f ' reduces to 

T' = 2if' f)r-<r T)r’-<r r')T 09 ) 

Note that if f = (f, both a and t are constant, otherwise, they vary with s. Having come 
this far, it is possible to define a very good (but not yet elegant) method of solving the ori- 
ginal set of differential equations. Using s as independent variable, and realizing that X and 
V can be arbitrarily chosen to be 0, integrate the set 

t’ = r 

r" = aT+f-t 

r" = ar +(L+f $ ( 20 ) 

o' = 2 (r-f) 

t ' = 2(f ' • f y- (r ■ T y - <r • r ’ )? 

The redundant equations are integrated because it has been shown that it is more accurate 
to integrate them than to calculate r, a and t from the integrated values of T and f , In 
fact, as a general role, which shall be applied again later, if at any level of knowledge (i.e.. 
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at any level of integration) integrated variables ace combined in a single operation, it is 
more accurate to integrate a variable defined to be that combination. Alternatively, opera- 
tions involving identities at a given level are probably equally accurate, except that an iden- 
tity with the fewest number of integrated variables will most likely be the most aocmate. 
The fundamental problem with integrating eqs(20) is that * is being integrated directly, and 
no tee is being made of the fact that a conic solution could be used to remove the effect of 
the central force term. Ideally then, one would compute a conic solution, and then compute 
the variation in the constants of that solution caused by the perturbing acceleration, f. This 
approach will yield greater accuracy than the direct approach, and is known as variation of 
parameters. Different variation of parameters methods vary different parameters, but they 
all begin with a solution to the unperturbed problem. A very elegant conic solution can be 
developed by expanding t and ? in Taylor series about the initial point, and assuming that 
f = Q = f = ft. The first few of the higher derivatives needed in the expansion for t are 
given below. Note that a pattern quickly develops, so that any of the derivatives could be 
written down by inspectio n 

t 9 — r 0 


to = r. = or. + p 
to = or. 

t 0 '”" = ar 0 " = o 2 r 0 +a*i 


( 21 ) 


etc. 


The Taylor series expansion for t is, of course 

t = t 0 + to's + + *o'" 4 r + 


2 ! 


3! 


Using eqs(21) in eq(22) and collecting coefficients results in 

s 3 jr 5 3 s 7 . 

t=to + r a (s + a— + a 2 — + a 3 — + . . . ) 

+ r.'(| + oil + a 2 |+... ) 

c 3 t 7 

+ M^ + «7r + o2 77 + --- ) 


3! 5! 


( 22 ) 


(23) 


Following the same approach, 

K" = aK -t 

n "" = ar o " * a¥. - <?*) 

K""' = o¥ 0 ' 

r 0 """ = a¥/' = - a* 

etc. 


Note how the linearization resulting from the introduction of the Laplace vector makes the 
computation of these derivatives trivial. The Taylor series expansion for T is, of course 




(25) 
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Using eqs(24) in eq(25) and collecting similar coefficients results in 

r-niwad. ■*£♦«»£♦... ) 

The results obtained thus far can be simplified greatly by defining the following transcen- 
dental functions given by Goodyear [7] which are a farm of Stumpff functions [8]. 

S„ * 1 + o4- + a 2 — + a 3 — + . . . 

° 2! 4! 6! 

5 1 = J + a^ + a2 ^- + a 3 ~ + ... (27) 

S 2 S* t S 6 

Sl *2T + °« + ' 

p r 5 s 5 it 7 

Sj E 3T + a sT + 7f + • • • 

A number of identities can be established by inspection erf the S series. For example, 

S 0 = l + aS 2 , 5 ! = s + aS 3 , etc. (28) 

Also, the derivatives with respect to s can be deduced by inspection (recall that when f is 
zero, a is constant ) 

S/ = oS, 

Si' = 5* (29) 

Sj =Si 
S 3 ' = S 2 

A summary of the identities used in this paper are given in Appendix 1. 


Using the new transcendental functions, the previous results can be written 
t =l 0 + r a Si +r/S 2 +\lSi 
r = /' = r a S 0 + r„'Si + pS 2 

r' = Qi+or 0 )S| + r o ’S 0 (30) 

? = +K'St -ts 2 
r’=(af 0 -ttfi + r 0 's 0 


This then is the solution to the unperturbed conic case. Eqs(30) may be used to determine t, 
r, r', f*, and T’ as a function of s and the initial conditions, t„, r 0 , r 0 \ ?„ , T 0 a, and 
These equations are valid on any conic. (i.e., elliptic, parabolic, or hyperbolic), without 
exception and without singularity. Also, very efficient routines exist for computing the 
Stumpff functions for any value of s or o. A very elegant method of computing the series 
is outlined in [7]. 

The appearance of the terms p + ot„ and a T„ - £ in eqs(30) result from taking the deriva- 
tive of S 0 with respect to s. In the variation of parameters method that is to follow, all 
parameters that are constant now, for the unperturbed case, will be considered to be vari- 
ables. This means that the sum and differences of integrated variables will be used to 
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compute Y ' and r\ This is not a particularly good thing to do for the following reason. 

Each of the variables, t, a, and Y„ will have some integration error associated with them. 
Forming the difference as indicated in eqs(30) could lead to a magnification of this error 
due to a loss of leading significant digits. This problem can be avoided by defining 
» aY c - 1. and y„ ■ 4 + ar 0 . Then, with the use of eqs(28), eqs(30) can be re-written 


as follows 



t = i 0 

+ r„Si 

+ r„ # S 2 + pS, 

r = 1 ’ = 

r,S, 

+ r/Si + |xS 2 

r' = 


r,'S 0 + y e Si 

K 

+ K'S 

1 + ^0^2 

Y’ = 

K'S C 

+S c Si 


This form of eqs(30) should be as accurate as possible since each member depends on the 
fewest number of (what will be) integrated parameters. 

Variation of Parameters 


The essence of the variation of parameters method is to develop equations for the variation 
of the parameters used in the un pert ur bed solution. In the perturbed case, when f, $, and 
7* are not zero, eqs(30) or eqs(30a) can still be used to compute t, r, r', Y, and Y ', provided 
the constants or parameters of the motion, Le^ i # . r a , r e ’,Y 0 ,Y e \ ft, , y 0 , a, and P arc con- 
sidered to be functions of the independent variable, s. Because the parameters should vary 
so much more slowly than the state, a significant improvement in the accuracy of comput- 
ing the state can be expected. 

Differential equations for these parameters are developed below. It doesn’t really matter 
whether eqs(30) or eqs(30a) are used in the development, the end results differ only by 
identities. Therefore. eqs(30) will be used now, and the computation of £ and y will be 
postponed until later. 

Defining ft * Y a and ft = Y 0 \ eqs(30) can be used to compute Y and Y ' as 

Y = +ftSj -fS 2 (31) 

and 

P'Map-ttfi + ftft (32) 

provided the derivative of eq(31) with respect to s is the same as eq(32), and the derivative 
of eq(32) with respect to s is the same as eq(13) which is repeated below 

Y" = aY+f-t (13) 

The imposition of these two criteria are sufficient to develop the variation of parameter 
equations for p and ft. Differentiating eq(31) 

Y' = p% +ft'S, -t'S 2 + J*S 0 ' + ftS,'-« 2 ' (33) 

In this case, the computation of the 5/ variables must take into account the fact that a is 
no longer constant. Denoting 4^ by ( )° and recalling eqs(29) 

V --- a$, + 5 “o' 

Si - S c + S]°o' (34) 

5 2 ' = S, + S 2 a <x' 

Sj' = S 2 + S“a' 
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The evaluation of the partials with respect to a will be addressed later. For now, using 
eqs(34) in eq(33), and requiring that eq(32) still be valid, results in 

ff = p 'S 0 + 'S, -?’S 2 + [jJ5“ + ffs,° - tS“jo' (35) 

Similarly, taking the derivative of eq(32) with respect to s, and requiring that eq(13)still be 
valid, and using eq(28), results in 

f = (o}f'-t')Si + F'S. + [(aJf-Wr+Ks^ + pS.ja' (36) 

Now, multiplying eq(35) by S„ , and eq(36) by Si and subtracting, results in 

^tf=(5 < , 2 -oSj 2 V / -(S 2 S 0 -S I 2 )r'+ (37) 

[(S C S“ - otSiS® - S?)p+ (S 0 S,° - S,S“)?- (S„S? - S,S,°)^a' 

Using the identity 


Sg-aSi = 1 


(38) 


which is proven in Goodyear [7] by comparing terms in the expansion, and eq(28), it is 
easy to show that 


S i — S C S 2 = S 2 


(39) 


Also, using the identities 2S “ = sS t and 2S“ = sS 2 -S 3 , it is easy to show that 
S 0 S 0 a - aSiSf = Sf/2. Goodyear [7] proves the identities S,S“ - S 0 S “ = fsiS 2 + S 3 ]/2 
and S]S] 0 - S a S 2 = S 2 H. Using these identities in eq(37) results in 

p ' = -S,f - S^T ' + [s?p + (S,S 2 + S 3 )$ - Sir] y (40) 


Now, multiplying eq(35) by oS,. and eq(36) by S 0 and subtracting, results in 

S 6 f= (S 2 - <*S , 2 )ff ' - (S e - oS^S ,r * + (41) 

[(S,So + ct(SoS, 0 - S 0 “S i))JJ + (S 0 “So - aS, a S,)|r - (S,“S 0 - aSfS^a' 

Using the identities given previously, and the following two. also given by Goodyear [7], 
2 S“ + S 0 S 3 = SiS 2 ,and2S 2 + SiS 3 = S 2 , it is easy to show that 

H ' = So f + s,r ' - [(s + SoS,)jj + s , 2 ? - (S,S 2 - S 3 )r] f <«) 

Following a similar procedure, define p * r a and P = r 0 \ Then r and r' from eq(30) can 
be written 

r =pS e + PS, +uS 2 (43) 


r' = (n + ap)S, + pSo (44) 

Taking a derivative of eq(43) with respect to s, and insisting that eq(44) be valid leads to 

SoP' + S,p' + (pSo“ + PS? + pS 2 °)a' = 0 (45) 

Now, taking a derivative of eq(44) with respect to s, and insisting that eq(15) and eq(43) be 
valid leads to 

oSip' + SoP' + [(p + oq>)S“ + PSo" + S,p]a' = r@ (46) 

Now, multiply eq(45) by S 0 and eq(46) by Sj and subtract to get 
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(5 2 -aS,V + 

[(5.5“ - S , 2 - o5,5,“)p + (5.5, “ - 5,5“)P + (S 0 $2 - S,S{*)n] o' = - * • 05, (47) 

Using the identities given earlier, eq(47) can be reduced to 

P' = - f 05, + [s?p + (5,5 2 + Sj)P + p5|] y (48) 

Now, multiply eq(4S) by aS , and eq(46) by S. and subtract to get 

P'(So - oS?) + (49) 

[(S,S. + 0 ( 5 , 5 “ - 5|5,“))p + (5.5,® - 05,5,“# + (5.5, “ - o5,5 2 “)p)o' = f$S 0 

Using the identities given earlier, eq(49) can be reduced to 

P' = t 0S. - [(r + 5.5, )p + 5?P + (S,S 2 - 5,#] y (50) 

The last step is to develop a differential equation for t, . Defining x * /, and taking a 
derivative of t in eq(30), and still requiring t ' = r = pS, + PS, + pS 2 leads to 

x' + 5,p' + S#' + (S,“p + S“P + 5“p)a' = 0 (51) 

Now, substituting eq(48) and eq(50) into eq(Sl) and simplifying leads to 

t' = f 0 S 2 - [(5,S 2 - S 3 * + 5|p + (S,Si + sS 4 - 35j)pj ^ (52) 


Finally then, instead of integrating the set, eq(20), the set 

o' = 2 ur’f) 

f=2trTy-<rTy‘~<rr')f 

P ' = -5,f - Sy ' + (s?|f + (5,5 2 + 5j)? - s 2 2 *| 

' = S.f + 5,t ' - [(j + 5.5,# + 5?fl - (S,S 2 - Sjje] y (53) 

p' = - f ■ 05, + [s, 2 p + (S,S 2 + Sa)P + MS 2 2 ) y 
P f = f [(r + SoS,)p + 5?P + (S,S 2 - S 3 )p] y 

X' = f 0 S 2 - [(S,S 2 - Sj)p + S 2 2 p + (SjSj + j5 4 - 35 s )u] y 

could be integrated, with, of course, eqs(30) used to compute /, T*, Y ', r, r\ 

As was mentioned earlier, eqs(30a) are probably a better set to use to compute 
l,T,T',r , r ', since combinations of integrated variables are not used together when multi- 
plying a single Stumpff function. If this is done, % and y must be integrated directly. Since 
there is no need to integrate both t ' and $ ', t might as well be removed from eqs(33). 
This can be done easily by first examining only terms involving p and t in p '. 

S 2 |J- S 2 (ap-$) = (S, 2 - aS 2 # + S 2 8 

But from Appendix I, 5, 2 - a S 2 = 25 2 , hence p ' becomes 

p' = -Sif-SiT+ [ 252 ? + (S,S 2 + S 3 )P + 5 2 2 S] y 

Tuming now to the coefficients of p and P in ft ' leads to 
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(5 + S 0 S j)jt (S jSj ■“ •?) s (j + 5 0 S] ** aS x S 2 + oS^lft + (S jS 2 — 

Using cqs(28) to replace s and 5 a , leads easily to 

P ' « S 0 f + S x t ' - [2S,Jf +.*, 2 fl + (S iS 2 - Sjtfj ^ 

Since S s ap - P, it’s derivative is simply 

$ ' = op + ptt' - 1 ' = - aS , f - S, * ' + [2S 0 p + o(.9 ,S 2 + S 3 )ff + aS 2 2 J>] y 

Considerable numerical experience seems to indicate that, contrary to what one might 
expect, computing y makes no appreciable contribution to accuracy, consequently it will not 
be considered further. 

In summary then, in lieu of integrating eqs(53), one could integrate eqs(53a) below 
a' = 2( t'T) 

p ' = -S,f - SjP ' + [ 2 Sjp + (S t S 2 + Sstf + S 2 2 $] y 
S ' = ap ' + pa' - - aS , f - S,P '+ ( 2 S, p + a(S jS 2 + S 3 # + aS f$] y 

$ ' = S 0 f + S,t ' - [ 2 S jp + S, 2 £ + (S,S 2 - S 3 j5| ^ (53a) 

p' = -rgS 1+ [s, 2 p + (S,S 2 + Sa)P + M^l] y 
P' = f • & - ((J + S 0 Si)P + S?P + <S,S 2 - 5 3 )u] y 
** = r £ S 2 - [(S,S 2 - Sj)p + S 2 2 p + (SjSj + sS 4 - 3S s )p] ^ 

The integration of eqs(53a) where t ' is computed but not integrated, allows the use of 
eq(30a) for computation of ihe state without forming f from other integrated terms. Also, 
it is probably equally accurate to integrate either form of $ since they differ only by iden- 
tities at the same level of knowledge. 

The Optimal Use of Fictitious Time 


The interesting point to note in all this is that, whichever set of parameters is used, 
f*. 7", r , , r ' and t are all computed from the solution set p, }J, p, p. a, and t (or and t, 
using the current value of s. 

Most people seeing these complex equations for the first time tend to be overwhelmed by 
the derivation and don’t really see through all the algebra just what is happening physically. 
However, the equations developed above can be interpreted schematically so that an intui- 
tive feel for wha: variation of parameters is up to can be obtained. 

A real understanding of variation of parameters may be obtained by assuming that a "true", 
i.e. perturbed, trajectory path is generated by integrating eqs(20) forward from t 0 . Note that 
this "true" path includes all perturbations and is a function of s. At every point along this 
"true" path, T and t ' are known, and may be used to define a two-body conic, i.e. the 
osculating conic. Now, assume that at each point along the "true" path the osculating conic 
is used to coast backwards a distance -s, and that having arrived at -s, the osculating 
parameters, called "of, are re -computed such that now a forward coast of +s brings us back 
to the "true" T* and T '. One could imagine then the unfolding of "Ct" as s u increased from 
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zero, and that given s, one would know *tr, and b- could compute the "true' r and t \ 
But this is exactly what eqs(30) or eqs(3Qt) imply. Th equations for a and t ae truly 
osculating, and they are the same on the ‘true’ path and after masting back -s. The vari- 
ables p. ff, p. and P may be thought of as functions of -s, since they are used to coast for- 
ward +s to define Y, Y \ r, and r' along die ‘true" path. It should be obvious that there are 
an infinite number of these that could be used depending on the choke of -s. The simplest 
choke is to choose s = 0. which implies 

p = r *m, ? p' = r M * 



ft ! 

ft ,f 


r = r'=pS, +p5, + |lS 2 
r’ - (p-*-ap)5 ; + 

r=ps c +$Si-ts 2 
F*' = (ap-m, + fc 0 


r 


Figure 1. Schematic Interpretation of Variation of Parameters 


These ideas are depicted in Figure 1, in which the "mie" trajectory unfolds to the right, 
while the variation of parameters integration occurs in the "IT space. The problem with 
keeping the initial point fixed can now be seen rather easily. If the perturbations are 
periodic, there will not be much gross motion in the "l t" space since the basic orbit is not 
being changed excessively. However, when persistant perturbations occur, such as long 
term drag or low thrust, there will be continual motion in the "0 T space. If the "true" path 
changes such that a wildly different conk is needed to coast to that point from the initiol- 
pouu , rapid escalation in the "(T space is possible, and trea^ .own of the method can 
occur. In two specific cases, low thrust to escape and re-entry to the ground, breakdown of 
the original K-S, BGI4, and BGI45 methods occurred. Now, consider what happens if the 
‘0T space is not hooked to the initial point, but rather moves with the trajectory. Now the 
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only motion that occurs in "(T space occurs during a single step, and them is no accumula- 
tion of effect that must be coasted out with an accumulated independent variable. This is 
also shown in Figure 1 . Consequently, by sliding the origin, the escalation can be reduced 
considerably, or eliminated altogether. Another advantage of not having s increase mono- 
tonically, is that previously, aside from the initial point, nothing could be done to improve 
the right hand side of the variational equations This point has not been addressed yet, 
although it was alluded to earlier with the introduction of X and St. The idea is this. If the 
variational equations change slowly then larger integration step’ may be taken. By choosing 
X and V carefully it is possible to drive the forcing terms on the right tend side of eqs(53) 
to zero, at least at t a . However, if s is set to zero at the beginning of each integration step, 
X and ? may be chosen again and again, so that the right hand side can constantly be set 
back to zero. Also, if s - 0. all the 5, arc zero except 5, = 1. This means that 
t, jj, jj, p, P cm be reset io exactly t , f\ T '. r . and r\ with absolutely no error introduced. 


Note trc. this is not the same as a re-initialization, since a and t are not altered (other than 
changing X and V ). ^-initialization would result in a loss of accuracy which can be 
observed quite readily. In fact, the difference in the energy and radius as integrated, com- 
pared to those computed from t and t ' is used as a measure of integration accuracy. Con- 
sequently, advancing the origin. i.e„ seeing s - 0, accomplishes several objectives. It 
makes the program faster. It removes sensitivity to a distant initial point It removes the 
need to worry about secular terms, since the independent variable is always near zero. It 
allows the free parameters X and V to be chosen at earn step ir order to improve, accuracy 
and increase step size. Since the independent variable, s, is going to be to r-ro at the 
beginning of each step, it is interesting to look at the differentia! equations for the variation 
in the parameters, and see if a choice can be mad; for X and that will make the right 
hand side zero. The hope is that this will cause the parameters to change more slowly, and 
consequently allow larger integration steps. Although only numerical experience will indi- 
cate what strategy, if any, is beneficial, what follows assumes that minimizing the right 
hand side of eqs(53) is the best strategy. Examining eqs(55' shows that P' will be zero if 
r Q is zero. Choosing X to cause tlus to happen leads to (recall that f - Ft ) 

X = r-f (55.1 

Since a* is continuous even if X changes, from cqs(12i. u follows that 

Q* = Gtfie ~ *m - f**™- ~ X^. (.56 ! 


Consequently, 


ou=«Ui + ex*, - Ku) ' 57 > 

Again examining eqs(53) shows ft ' will be zero 3t s = 0 if T is zero. Consequently, 
choose 

V^w^-r 2 /* (58^ 

Also, since the perigee vector, J?, is continuous even if V charges, from east 12). it tollows 
that 


W5 — ^cld ^old ^ 


(59) 


Consequently, 


- t'oiz 


Void) 


(60; 


If and y are used, and p and p are updated, it is w'orse. since K must be jumped to keep T 
constant, and y must be jumped to keep p constant. If thereafter. X is changed, .hey will 
both be jumped again because a will change. If further, V is changed. S must be jumped 
an equivalent amount. This process can be thought of as jumping from surtace to surface in. 
multi-dimensional space in order io keep die parameters changing as slowly as possible. 
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Note that X^ and are held constant for the duration of an integration step, and are 
only changed at the beginning of the next integration step (when they become 'old'). Nate 
also that neither X nor V must be changed, it is just possible to change them if numerical 
experiments indicate that there is some benefit in doing so. 

Relation to Bond and Gottlieb’s Fourteen Element Method - BG14 

The objective of this paper has been to present a simple derivation of a variation of param- 
eters method, using tarn *sal variables and fictitious time, and to point out that the optimal 
use of ficuaous tune is to set it to zero at the start of each integration step. In other words, 
advance the origin. Advantages include faster running time and greater accuracy, secular 
terms are not important, and there is no build up of sensitivity to a distant starting point. 
Also, it possible to choose constants to drive the right hand side of the variational equation; 
to zero at the beginning of each integration step. 

The introduction of the V associated with t is also new in this paper. The use of the 
Keptcrian energy constant, a*, in this paper was motivated more by simplicity of presenta- 
tion than any desire to use the "best" energy constant. An observation of eq$(53) or 
eqs(53a) shows that a' is a major contributor to the right hand side of those equations. 
Consequently it seems reasonable to try to use for a a constant of the motion that, in bet, 
remains constant for the greatest possible range of pertisbing accelerations. Total energy 
would be a better choice than Keplerian energy, since total energy is constant in the pres- 
ence of a perturbing acceleration that is due to a time invarient potential. The Keplerian 
energy was used in [3], and the total energy' was used in [9], However, the Jacobian 
integral of [4] is best of all, since it remains constant in the presence of a perturbing 
acceleration resulting from a time varying potential. Therefore, in lieu of defining a* as ir* 
eo(9), define die Jacobian integral 

a « r ■ r- & + 2(V(r*i) - 0 ■ ?) (61) 

'r 

where t = t x t is die angular momentum, and 0 is a constant vector which characterizes 
the time rale of change of the potential. Now, defining asffl?, solve for the Keplerian 
energy as 

r - r- = a, - 2 (V(?V ) - a) (62) 

r 1 

and substitute into eq(8) to get 

T“ - r 2 f - a, - 2(V{f'j) - a) jr* (63) 

as before, add and subtract IT and V to get 

r " = a/+ r 2 f - 2(V(r.t) -a? -pjt + kr-xr + v-? (64) 

Now define 
a ss a.j + X 

?E)Jjf + V ( same as before ) (65) 

(T *r 2 f -xr-2(V(rj)-oy 
f*Q+? 

With these definitions, eqs(53) and eqs(53a), given previously, are still valid with the 
exception that now, from eq(61), a' becomes 

r > 

o' = ra, = r\2? - ? + Qf + 2V' > * r - 20 • t* 2V' : j (6b) 

1 i 
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where superscripts an V denote partial derivatives by r and time. Recall from eq(5) that 
r - t ■ F7r. Also, [4] derives the fact (hat V‘ = - 0 ■ f x V'’. Using these, it is easy to 
show that eq(66) reduces to 

a' = 2(f‘ , -raxO? (67) 


Note that now a' depends only an $ and not / . Also, because a was introduced into the 
problem, it is integrated from it’s initial point rather than being calculated. This is done 
because it is more accurate to do so. The differential equation for a is 


<f =rd=rJ3rxr*=r5Jr‘x/* 

(68) 

The only other difference is that now the optimal X is 


\ = r-T +2io-v(fj)) 

(69) 

and Vf is 


?=-£? = (X + 2(V-a))r*-r 2 / 

(70) 


The major advantage in changing energy constants in the derivation is that greater accuracy 
can be attained by using an energy integral that is constant for a greater variety of perturb- 
ing accelerations. The cost for this is the introduction of the differential equation for a. 
However, the advantage that accrue? from the use of the Jacobian integral is more than 
worth the extra differential equation. 

Numerical Results 

The formulae developed above have been coded in the Ada Simulation Development Sys- 
tem (ASDS) {13], and tested on a number of problems. One problem run was the 288 day 
J2 plus moon problem given on page 121 of [10]. This problem was run with and without 
and advancing the origin. The answer in [10] was obtained in both cases, when SO steps 
per rev were used in conjunction with Runge_Kuua_Fehlberg_78 fixed step. When 30 steps 
per rev was used, advancing the origin and choosing V resulted in a slightly smaller root 
mean square error. Although this problem contains a perturbing force outside the potential 
term, it is not a particularly severe test of the algorithm. Two severe test cases were run. 
The first was nominally a reentry to the ground. It was actually set up as a fixed final time 
problem, with the final time being chosen to bring the spacecraft essentially to the ground. 
Essentially is used because the final time was truncated to eight digits, 29.422577 days. 
The reentry problem: The reentry problem starts at t = (6677833.0 -62810.0 -2730!.0Xm), 
t = (79.0 6821.0 3627.0Xm/s). The ballistic number was 78.00cg/m**2). The gravity 
model was GEM10, 2x2. The atmosphere routine was Jacchia ’70, with the density held 
constant below 90km. The results of this problem are given in Table I. All runs were made 
with Runge Kutta Fehlberg 7/8 fixed step, and X and V were held constant except in those 
runs where X or v appear. The second severe test was a low thrust spiral from a circular to 
a hyperbolic orbit. The low thrust problem: Steve Sponaugle {111 first attempted to solve 
this problem using the K-S method [2], He noted that the program could not start in a cir- 
cular orbit and integrate all the way to hyperbolic. However, he also noted that if he re- 
initialized about two-thirds of the way out the program would integrate on through the para- 
bolic region, albeit poorly. It was this one clue coupled with the desire iO initialize prop- 
erly, that led to the sliding origin concept. The low-thrust problem starts at t = 
(6783140.0, 0.0, O.OXm), t - (0.0, 6736.774, 3657.770Xm/s). The thrust and flow rate are 
122.36657(N) and 2.4955962e-03(kg/s) respectively. The initial mass is 375000.0(kg). 
Thrust is along the velocity vector. The gravity model was GEM-9. J2 only. Several runs 
of 240 days were made. The first two were made with the sliding origin method of this 
report using 80 steps/rev. Final r* = ( -634757780.495, 898972406.577, 20257150.2835)(m). 
This is shown in Table II. The other runs in the table compare the BGl-ib and BG14e slid- 
ing origin (sf) methods to each other as well as to the BG145 method in [12]. Note that 
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again the simple removal of secular terms failed to cure the instability. The addition of the 
sliding origin to the BG145 of [12] results in a method whose rharw-remtirs are identical 
to those of the BG14S of this report. Further, the effect of causing X to change at each step 
is shown to have a dramatic effect in the low thrust problem, but very little for the reentry 
problem. Also, changing V at each step as well as X actually caused things to get slightly 
worse. Note that in every case the c farm is more accurate than the 6 farm, although not 
by much. 

Conclusions 

An examination of Tables I and II indicates that simply removing secular terms does not 
correct the instability in the original BG14. However, when the sliding origin was a AW to 
Bond's BG145 method (BG14SQ in the tables), that formulation and the BG14&Q formula- 
tion herein give essentially identical results. The BG14cQ method gave essentially identi- 
cal results for the low thrust problem, and, contrary to what was expected, gave slightly 
superior results for both the reentry pobiem and the low thrust problem. Also, choosing X 
at each step improved the low thrust problem, but had little effect an the reentry problem. 
Previous studies had indicated that the 5Q form of BG14 was slightly superior to the eQ 
form. The results here do not show that. However, previous studies used 4th order 
Runge Kutta rather than the more accurate 8th order used herein. The conclusion one 
should draw is that the methods are essentially the same and perhaps other criteria such as 
simplicity and speed should enter into the decision for a choice of formulation. Since the 
methods differ b> only a few terms, one could code them both or arbitrarily choose one. It 
is suggested that studies be carried out with specific types of problems to see if X or V 
should be changed at each step, and whether the SQ or the eQ form gives better results. It 
is obvicus from [12] that either the change of variables or the removal of secular terms 
improves the original BG14 [4] if the origin remains fixed. However, since that method 
failed to solve either the low-thrust problem or re-entry to the ground, one must conclude 
that the real culprit is the fixed origin and not secular terms per se. When the sliding origin 
was added to Bond’s S method, the instability vanished, and it solved the problem getting 
the same answer as obtained by the formulation represented by eqs(53a). The formulation 
represented by eqs(S3a) is essentially the same as in[12], with the exception that y does not 
appear on the right hand side of the differential equations and no particular effort was made 
to remove secular terms. Since eqs(53a) and the formulation in [12] (with sliding origin 
added) give essentially identical results, one must conclude that it is not the form of the 
right hand side of the differential equations that really matters, it is rather the choice of 
dependent variable that actually makes die difference. Finally then, the optimal use of ficti- 
tious time, i.e., moving the origin renders secular terms innocuous, since they can never 
grow, and also overcomes the sensitivity to the initial point that results from long term per- 
sistant forces. Any element method using fictitious time can be stabilized by setting ficti- 
tious time back to zero as indicated in this note. The addition of this concept to any of the 
Jacobian <ntegral methods, ix., BG14 in any of its forms ([4], [12], or as developed herein), 
makes them highly accurate, long term propagators. 
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Appendix I 

Useful identities, which are either given in [7] or which can be derived from them are 

S, 3 -aS? = 1 
S 0 = 1 + aS 2 
S; = s + aS 3 
Si - S C S 2 - S 2 
S 1 = Z S 2 + aS 2 2 
2 S? = sS x 
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2ST =iS 2 -S 3 = S,S 2 -S 0 5 3 
2S° =sSj-2S 4 = Si ~s t s 3 
2Sj « sS 4 - 3S s 
2 (S t Sf-S 4 Sn = S l S 2 + S 3 
2CW - SyS?) = S^J 
2(W-$ # 5 2 a ) = J 2 2 


Classical Stumpff Functions 


If one defines the variable z ■ cu\ and forms the following transcendental fimcfions 

.2 ,3 

■+■ - 

4! 

* 

z* 


1 z z* Z J 

c '*7jr + 2T + 7r + 7r + 


6 ! 

.3 


• I 


1! + 3! + 5! + 7! + ' 
1 z z 2 z 3 
2! + 4! + 6! + 8! + ' 


1 z z* 

C3 *^' + sr + 7T + or + - 


d ii 

3! ' 5! ' 7! + 9! 

It is rather easy to show that the c’s and S’s are related by 

S„ - c„ 


5 1 = sc i 

5 2 = s 2 c 2 

$3 = ^3 

5 4 ~ s*c 4 

5 5 = s s c s 

The c’s are known as Stumpff functions and are used in (4] except that a in [4] is the nega- 
tive of a in this report The S form of the Stumpff functions is used by Goodyear [7], and 
tend to make the equations seem less cluttered. 


APPENDIX II Multiple Vehicles and Constant Time Steps 

Since the independent variable is fictitious time, and since fictitious time is defined by 

— = r, it is obvious that smeejn general, the radius of one vehicle will not be identical to 
ds 

the radius of another vehicle, time will unfold differently along the two trajectories. This is 
inconvenient when looking at relative motion between two vehicles, since almost all 
relevent measurements require time simultaneity. It is possible to gain time simultaneity 
the following way. First of all, one of the vehicles must be chosen to be the master vehi- 
cle. This vehicle will have fictitious time, s, truly independent. Assume this is vehicle one. 
We would like time on the second vehicle to unfold at the same rate it does on the first 
vehicle. That is to say, we would like 

dt - r,d j! = r^ls 2 (H*l) 
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or 


ds 2 - — -dsi 
r 2 

This can also be written 

ds 2 >-i 
ds i r 2 


GI-2) 


ai-3) 


Therefore, if eqs(53) for the second vehicle are multiplied by — , the independent variab.- 

is shifted to s\. In addition however, whenever s 2 is needed in order to compute r, r\ etc., 
it must be computed as the integral of eq(II-2), i.e., as 

*2 = J — ds (n-4) 

J r 2 


Now that this is understood, it is easy to see that vehicle one is not even reeded if a con- 
stant time step is desired. In this case pretend that vehicle one exists, and that this ima- 
ginary vehicle is in a perfectly circular orbit, and that it is being integrated with a step size, 
ds, that is chosen to be 

ds = — (II-5) 

n 

where dt is the desired time step, and r, is just a reasonable value for the radius of the ima- 

r, 

ginary orbit. Then, eqs(53) can be modified by multiplying by — , and j] is computed 
from 


*1 = f — ds (n-6) 

j ri 

In this way a single vehicle or many vehicles can be integrated at a constant time step size, 
with no iteration what-so-ever and with great accuracy. Note that the equation for t is no 
longer needed since time is known absolutely. If t is integrated time can be computed and 
used as a measure of integration accuracy. 

A test case for a near circular orbit perturbed by a 4x4 gravity model was integrated for 30 
days using both the standard BG14 and this constant time step approach. The step size was 
varied from 25 steps per rev to 8 steps per rev. In every case, the rss error away from a 
high precision run was slighdv less with the constant time step method. An interpretation 
of this result is that the error made in integrating for s results in less error than that result- 
ing from the computation of time using X. The method above is convenient and seems to 
be quite accurate, at least for near circular orbits. However, since the computation of ficti- 
tious time involves the integral of 1/r, this approach may not be as accurate for elliptic 
orbits as the standard method that computes time analytically. A variable step size 
approach would be required as a minimum. In addition, in the absence of perturbations, an 
error would be introduced by the integration of 1/r, that would not be present with the stan- 
dard method. A variation of this first approach would be to integrate T, and solve for s 
from the time equation at every function call rather than integrate 1/r. This would add 
iterations for s, but not in terms of function calls. The advantage of this approach is that 
no error would be introduced by the integration of 1/r. Information from the next section 
can be used to speed convergence. 

A SECOND APPROACH TO CONSTANT TIME STEP 
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A second approach to stopping on a specific time (with s as independent variable) can be 
developed by assuming that s is a function of t, and forming the Taylor sen.'S 



s=4 + i* + ^$" I+ ¥$ 4, ’ + - 

ai-7) 

The derivatives of s that are needed come from the fact that s is defind by 



dt 

ds ~ r 

(n-S) 

hence 

ds 1 

di r 

(n-9) 

and 

dh r 

dt 2 r 2 

(n-io) 

where 

r-r 

r = 

r 

(n-ii) 

also 

d 3 s r 2 r 2 

dr 3 r 2 + r 3 

(n-i 2 ) 

Now, 

r-r r-r r-r. 

r = + r 

r r r 2 

(0-13) 

Since the velocity vector, V = K and remembering the definition of r , r becomes 



v 2 t? r 2 

r = — + 

r r r 

(n-u) 

Now, from eq(l) 

r-f^r-f- f 

(D-15) 

Consequently, 

d 3 S 1 It 2 - 7\ 

— = — (3r 2 + ^ - v 2 - r ■/) 
dr 3 r 3 r 

(n-i6) 

Defining 6s — 

, there results 



As»S-|r* 2 + j(3r 2 +*- V J -r-/>8 J 

2 6 T 

cn-17) 


Whenever, 5 is a reasonable number, eq(II-17) can be used to estimate As to land on a 
specific time. Since a third order Taylor series is being used, the estimate of As should be 
very accurate, and the time should be hit with very few iterations, probably only one. 
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TABLE I - Comparison of BG14f£2, BG146Q and BG14S Methods 
Final Value of Position Vector - Reentry to the Ground (T = 29.422577Days) 


Method 

Gottlieb 

X (km) 

Y (km) 

Z vkm) 

BG14CQ- lOOspr 

2517.0799480 

-5737.8247524 

1188.3706594 

BG14&2 - lOOspr 

2517.0799488 

-5737.8247526 

1 188.3706592 

BG14tf2 - lOspr 

2523.8693524 

-5739.4783332 

11862762008 

BG14&Q- lOspr 

2523.9860855 

-5739.4725808 

1186.2269801 

BG145Q-8spr 

2524.905855 

-5739.6584337 

1185.919466 

BG14£Q-8spr 

2524.82146 

-5739.649530 

1185.949942 

BG14S£2-8sprX 

2524.90541 

-5739.658044 

1185.919490 

BG14fi2 -8sprX 

2524.82098 

-5739.64925 

1185.950022 

BG14SQ- 8spr X v 

2524.90195 

-5739.65710 

1185.920422 

Bond [12] 

Failed 

Failed 

Failed 

BG14S - lOOspr 

BG14&Q- lOOspr 

2517.0799489 

-5737.8247527 

1188.3706592 

BG145Q- lOspr 

2523.9860855 

-5739.4725807 

i 

j 

1186.2269801 ] 

! i 


RSS(km) 


7.41656 
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TABLE II - Comparison of BG14e£2, BG14&Q and BG 145 Methods 
Final Value of Position Vector - Low Thrust Spiral to Escape (T = 240.0Days) 


Method 

X (km) 

Y (km) 

Z (km) 

RSS (km) 

Gottlieb 

-634757.780495 

898972.406577 

20257.150284 


BG146Q - 30spr 

-634757.776479 

898972.409549 

" 

20257.147576 

0.00568 

BG14SQ-30jpr 

-634757.776507 

898972.409528 

20257.147595 

0.00564 

BG14£fl - lOspr 

-634690.244994 

899021.597745 

20211.855408 

95.0391 

BG145Q- lOspr 

-634690.245016 

899021.597729 

20211.855424 

95.0391 

BG14&Q- 8spr 

-634262.290049 

899331.417762 

19925.408943 

696.026 

BG14eQ • 8spr 

-634265.16938 

899326 609579 

19928.146539 

690.195 

BG14&Q- 8sprX 

-634788.200341 

898950.438100 

20277 J 18474 

42.695 

BG14e£ 2 - 8spr X 

-634788.196759 

898950.436222 

20277.51741 

42.693 

BG14&l-8sprX v 

-634788.648135 

898950.114106 

20277.818275 

43.324 

Bond [12] 




1 

i 

BG145 - 30spr 

Failed 

Failed 

Failed 

N/A 

BG14§£2-30spr 

-634757.776479 

898972.409549 

20257.147576 

0.00568 

BG14&G;- lOspr 

-634690.244994 

899021.597745 

20211.855409 

95.0392 

i 
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